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20  ABSTRACT  (Continued) 

the  complex  plane  when  a  parameter  of  the  equation  is  varied.  The  method  is  a  two-dimensional 
counterpart  to  the  one-dimensional  technique  whereby  the  change  of  sign  of  the  pertinent  function 
delimits  the  location  of  a  root.  The  complex  root  is  similarly  enclosed  in  a  nested  set  of  squares  of 
diminishing  size.  The  method  is  illustrated  by  a  typical  example,  the  dispersion  relation  for  the 
propagation  of  straight-crested  waves  in  a  homogeneous  plate.  Without  fluid  loading  the  propagation 
speed  is  real;  loading  the  plate  by  a  fluid  moves  this  real  root  into  the  complex  plane,  which  physically 
corresponds  to  the  appearance  of  radiation  into  the  fluid.  The  real  and  imaginary  parts  of  the  relative 
wave  speed  are  presented,  calculated  according  to  exact  elasticity  theory  and  thick-piate  theory,  for 
antisymmetric  and  symmetric  waves  separately  and  simultaneously,  as  a  function  of  the  dimensionless 
wave  number.  Flow  diagram,  source  program  listing,  and  computation  examples  of  the  FORTRAN 
program  are  given. , 
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COMPLEX  ROOT-FINDING  PROGRAM  WITH  APPLICATION  TO  THE  DISPERSION  RELATION 
OF  WAVES  PROPAGATING  IN  A  FLUID-LOADED  PLATE 


INTRODUCTION 

To  find  complex  roots  of  an  equation  one  often  starts  from  the  "princi¬ 
ple  of  the  argument",  which  can  be  stated  as  follows  [1]: 

"Let  f(z)  be  analytic  within  and  on  a  simple  closed  contour  C, 
except  possibly  for  a  finite  number  of  poles  inside  C.  If  f(z) 
does  not  vanish  on  C,  then 

(^r)  uc  arg  f(z)  =  Z  -  P  (1) 

where  arg  f(z)  is  the  total  change  in  the  argument  of  f(z) 
around  C,  going  in  the  counterclockwise  direction,  Z  is  the 
number  of  zeros,  and  P  is  the  number  of  poles  inside  C.  In 
counting  the  number  of  zeros  and  poles  a  zero  of  order  m  is 
counted  m  times,  and  a  pole  of  order  n  is  counted  n  times." 

If  it  is  known  that  there  are  no  poles  of  f(z)  inside  C  this  theorem  is 
obviously  useful  in  finding  areas  of  limited  extent  in  the  complex  plane 
where  zeros  are  located.  In  practice  problems  arise  in  choosing  the  step 
size  sufficiently  fine  so  that  no  zeros  are  missed  due  to  the  ambiguity  in 
the  argument  of  a  complex  variable.  This  is  discussed  in  reference  2.  To 
determine  the  location  of  the  zero  more  accurately,  one  might  go  through  a 
succession  of  smaller  contours — or  rather  use  some  other  approximation 
scheme,  as  for  instance  the  Newton-Raphson  method  [3]. 

In  this  study  a  different  method  of  root  finding  is  introduced.  This 
method  may  be  considered  as  a  complex  counterpart  to  a  familiar  root¬ 
finding  technique  for  real  functions.  In  the  real  case  one  looks  for 

intervals  where  the  value  of  the  function  changes  sign.  By  repeatedly  sub¬ 
dividing  the  interval  and  looking  for  changes  in  sign  of  the  function,  one 
may  determine  the  root  to  any  desired  degree  of  precision.  In  the  present 
routine  the  complex  root  is  approximated  by  a  sequence  of  squares  of  dimin¬ 
ishing  size,  within  which  the  root  is  known  to  be  located.  The  basic  algo¬ 
rithm  is  dependent  on  the  fact  that  in  the  neighborhood  of  a  zero  of  order 
one  a  complex  function  can  be  approximated  by 

f(z)  =  C(z  -  Z),  (2) 

where  C  is  a  complex  constant  and  Z  is  the  given  root.  To  obtain  a  first 
approximation  for  the  given  root  one  might  rely  on  the  principle  of  the 
argument,  Eq.  (1).  It  appears  more  attractive,  though,  to  make  use  of  the 
observation  that  in  many  practical  cases  the  function  f(z)  has  a  specific 
structure.  Often  the  function  f(z)  depends  on  a  parameter  r  in  such  a  way 
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chat  for  r  =  0  the  function  is  real  and  has  a  real  root,  which  can  be  found 
by  a  root  finder  for  real  roots.  An  example  is  the  dispersion  relation  for 
straight-crested  waves  in  plates.  Without  fluid  loading  the  propagation 
speed  is  real.  Fluid  loading  causes  an  imaginary  part  to  appear,  which 
corresponds  to  the  phenomenon  of  sound  radiation  into  the  fluid.  The  fluid 
density  features  as  the  parameter  r  in  this  case.  By  increasing  r  in  steps 
from  zero  to  the  nominal  valve  one  can  follow  the  march  of  the  root  from 
the  real  axis  into  the  complex  plane.  The  step  size  can  always  be  chosen 
small  enough  such  that  the  approximation,  Eq.  (2),  is  sufficiently  close. 

The  root  found  at  a  given  value  for  r  is  the  seed  for  the  determination  of 
the  root  at  the  next  larger  value  for  r.  It  is  clear  that  in  this  method 
one  cannot  find  those  roots  that  do  not  have  a  counterpart  in  the  real  case, 
as  is  true  for  some  waves  traveling  at  the  interface  of  a  plate  and  adjacent 
fluid.  On  the  other  hand,  it  has  the  advantage  that  a  given  complex  root  is 
identified  by  its  place  of  origin  on  the  real  axis.  Thus,  there  is  no  dif¬ 
ficulty  in  establishing  the  signature  of  a  root.  The  identification  of  com¬ 
plex  roots  from  a  random  set  in  the  plane  and  the  sorting  out  of  those  roots 
that  are  organically  related  can  pose  a  serious  problem  that  is  avoided  by 
the  present  method. 

In  the  second  part  of  this  report  a  typical  application  of  the  root¬ 
finding  method  is  described.  The  pertinent  equation  is  the  dispersion  re¬ 
lation  of  straight-crested  waves  propagating  in  a  fluid-loaded  plate.  Both 
the  dispersion  relations  following  from  exact  elasticity  theory  and  those 
that  are  based  on  thick-plate  theory  are  used. 

COMPLEX  ROOT  FINDER:  BASIC  ALGORITHM 

Let  F(z,r)  be  the  function,  the  root(s)  of  which  is  to  be  determined; 
r  is  a  parameter  such  that, for  r  =  0,F(z,0)  is  real  and  has  a  real  root. 
Suppose  that  one  has  two  points  and  zi  close  enough  to  the  root  Z  of 
the  function  F(z,r)  that  the  linear  approximation,  Eq.  (2),  is  applicable. 


X 

Fig.  1.  Relative  location  of  test  pair  z^.zt  ar.d  root  Z 


One  sees  in  Fig.  1  that  the  angle  0  =  0^  -  S2  between  the  vectors  from  Z 
to  the  points  z^  and  z?,  respectively,  in  the  complex  plane  is  in  the  first 
or  second  quadrant  when  the  loop  Z-'-Zj-*^  is  clockwise,  and  in  the  third  or 
fourth  quadrant  if  this  loop  is  counterclockwise.  The  angle  6  is  deter¬ 
mined  bv 

zj  -  Z 

-  =  '  r2  =  Tprr  •  (3) 

Since  the  sine  function  is  positive  in  the  first  and  second  quadrants  and 
negative  otherwise,  it  is  clear  that  the  location  of  Z  with  respect  to  the 
line  through  z^  and  Z2  is  determined  by  the  sign  of  the  imaginary  part  of 
(zA  -  Z)/(z2  -  Z) .  Another  way  of  interpreting  this  criterion  is  to  form 
the  vector  product  of  the  complex  numbers  (z^  -  Z)  and  (z2  -  Z) ,  considered 
as  two-dimensional  vectors,  given  by  (22  -  Z)x(zi  -  Z)y  -  (z2  -Z)y(zi  -  Z)x. 
The  algebraic  sign  of  this  vector  product  depends  on  the  orientation  of  the 
loop  Z->-zj->-z2,  and  also  has  the  same  sign  as  the  imaginary  part  of  (zj  -Z)/ 

(z2  -Z) .  If  zi  and  Z2  are  close  to  Z  the  following  approximation  will  be 
valid: 

zi  -  Z  C(Zl  -  Z) 

lm  -  =  Im  — — - — 

z2  -  Z  C(z2  -  Z) 

F(zl,r) 

~  Im  t'  /  — -  - 

F(z2,r) 

Im  F(zi,r)Re  F(z2»r)  -  Im  F(z2,r)Re  F(z]_,r) 

=  - 5 -  ,  (4) 

| F(z2  ,r) 1 “ 

where  Re  and  Im  indicate  the  real  and  imaeinary  parts  of  a  complex  number. 
Therefore,  the  algebraic  sign  of  the  expression  in  the  numerator  of  the  last 
part  of  Eq.  (4) 

S  =  Im  F(z^,r)Re  F(z2,r)  -  Im  F(z2,r)Re  F(z^,r)  (5) 

determines  the  location  of  the  root  Z  with  respect  to  the  test  pair  z^,z2. 
Using  the  words  horizontal  and  vertical  to  indicate  the  direction  of  the 
real  and  imaginary  axes,  respectively,  one  can  describe  the  core  of  the 
root-finding  procedure  as  follows.  From  an  initial  seed  Zg  close  to  the 
desired  root  a  vertical  test  pair  is  created,  namely  two  complex  numbers 
z\,  z2  with  the  same  real  part  as  the  seed  but  differing  by  a  small  step 
size  o0,  up  and  down  from  the  seed.  Thus, 

21  *  zs  +  iSQ 

and 

z2  “  zs  "  ifio- 

This  test  pair  is  moved  horizontally  by  steps  <$0  until  the  expression  S, 
indicated  here  by  Sjj»  of  Eq.  (5)  changes  sign.  A  similar  procedure  is  per¬ 
formed  with  a  horizontal  test  pair,  z2,z^,  determined  by 
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Z3  =  4(z1+z2)  -  6q 

and 

z4  =  1-s(z1+z2)  +  60  (7) 

where  is  Che  given  step  size  and  zj.z?  is  the  last  vertical  test  pair, 
after  SH  changed  sign.  The  pertinent  quantity  Sy  is  now  given  by 

Sy  =  Im  F(z3,r)Re  F(z4,r)  -  Im  F(z4,r)Re  F(z3,r).  (8) 

This  test  pair  is  moved  vertically  by  steps  of  the  same  size  until  Sv  changes 
sign.  Then  the  cycle  is  repeated  with  step  size  5^  =  5o/10,  etc.  until  the 
desired  precision  is  reached.  The  subsequent  vertical  test  pairs  are  formed 
by  setting 

zi  =  %(z3+Z4)  -  i5n 

and 

22  %(z3+Z4)  +  idn,  (9) 

where  in+i  =  6n/10,  and  z3,Z4  is  the  horizontal  test  pair  obtained  after  Sy 
changed  sign.  The  movement  of  the  test  pairs  according  to  Eqs.  (7)  and  (9) 
assures  that  they  are  at  all  times  as  close  to  the  root  Z  as  possible  in  view 
of  the  latest  approximation. 


The  initial  step  size  6a  is  determined  in  relation  to  the  change  in  the 
parameter  r  by  a  linear  approximation  of  the  function  F(z,r).  At  the  (i+l)th 
step  one  can  set  approximately 


F(zi+l,ri+i) 


(10) 


where  Ir  is  the  step  in  the  parameter  r,  and  this  linear  expression  is  used  to 
determine  5z.  The  partial  derivative  (9F/3r)^  is  approximated  by 


/  3F\  FCzi.r^Ar) 

V^/i" 

and  the  partial  derivative  (3F/9z)^  is  approximated  by 


(ID 


/ _  F(zj-h)z,rj) 


(12) 


The  quantity  iz  used  to  compute  this  difference  quotient  is  read  into  the 
program,  but  it  is  found  not  to  be  critical.  One  may  choose  it  to  be,  say, 
one  percent  of  the  nominal  value  of  jzj.  The  step  size  is  thus  given  by 


F(zi , r^+lr )Az 
F(zi+l1z,ri) 


(13) 


,t 
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This  step  iz  is  complex;  a  real  step  size  i0  is  found  by  taking  the  maximum 
of  the  real  and  imaginary  parts  of  Sz.  It  was  observed  in  practice  that 


Che  seep  size  found  in  this  way  is  not  always  adequate.  Therefore,  a  step- 
size  factor  is  read  into  the  program.  It  is  routinely  set  equal  to  one, 
but  if  the  program  does  not  converge  rapidly  one  may  enter  a  different  step 
size  to  remedy  this. 

APPLICATION  OF  COMPLEX  ROOT  FINDER  TO  THE  DISPERSION  RELATION  OF  STRAIGHT- 
CRESTED  WAVES  IN  A  FLUID-LOADED  PLATE 

Theory  of  Plane-Crested  Waves  in  a  Fluid-Loaded  Plate 

It  is  possible  to  find  solutions  to  the  wave  equations  for  infinite 
plates  that  are  straight-crested  waves  propagating  parallel  to  the  faces  of 
the  plate.  When  exact  linear  elasticity  theory  is  applied,  these  waves  are 

known  as  Lamb  waves.  In  this  report  the  development  of  the  theory  by 

Viktorov  [4]  is  followed.  According  to  the  outward  appearance  of  the  waves, 
one  distinguishes  antisymmetric  and  symmetric  waves,  which  corresponds  to 
different  parity  of  the  field  variables  as  a  function  of  the  coordinate 
perpendicular  to  the  faces  of  the  plate.  If  the  boundary  conditions  are 
the  same  at  both  faces  of  the  plate  (if  the  plate  is  loaded  by  the  same 

fluid  on  both  sides),  the  two  wave  types  may  occur  separately.  If  such  is  l‘ 

not  the  case,  the  two  types  occur  simultaneously.  . 

Since  exact  elasticity  theory  leads  to  considerable  complexity  for 
other  than  the  simplest  geometries,  approximate  plate  theories  have  been 
developed,  known  as  thick-plate  theories.  Thick-plate  theory  for  antisym¬ 
metric  waves  is  described  in  reference  5,  and  thick-plate  theory  for  symmet¬ 
ric  waves  is  described  in  reference  6. 

In  the  present  study,  the  root  finder  is  applied  to  dispersion  relations 
from  both  exact  theory  and  thick-plate  theory.  Moreover,  in  each  case,  one 
may  select  the  root  that  originates  in  an  antisymmetric  or  in  a  symmetric 

wave  for  an  unloaded  plate.  Another  option  is  to  choose  a  fluid  loaded  on  j 

one  side,  or  on  both  sides  with  the  same  fluid.  Although  strictly  speaking  ! 

the  boundary  conditions  are  not  satisfied  for  the  case  where  a  wave  of  either 

antisymmetric  or  symmetric  character  is  present  by  itself  in  a  plate  loaded  j  i 

on  one  side  only,  one  may  select  these  options  in  the  program. 

Dispersion  Relations  From  Exact  Elasticity  Theory 


The  waves  are  propagating  in  a  plate  of  thickness  2d.  The  symbol  c 
indicates  phase  speed,  k  is  the  wave  number,  and  the  subscripts  d  and  s 
refer  to  dilatational  waves  and  shear  waves,  respectively.  For  the  sake 
of  generality,  it  is  assumed  that  waves  in  the  fluid  are  incident  in  the 
plate.  The  waves  can  be  most  advantageously  expressed  in  terms  of  the 
potentials  <p  and  y  in  the  form 


and 


$  =  [As  cosh(qz)  +  Ba  sinh(qz)]  exp [ i (wt-kx) ' 
0  =  [Dg  sinh(sz)  +  Ca  cosh(sz)]  exp[ i(^t-kx) 


(14) 


where  x  and  z  are  coordinates  in  the  direction  of  and  perpendicular  to 
the  faces  of  the  plate.  The  phase  speed  c  equals  u/k,  where  ^  is  the  angu¬ 
lar  frequency,  k  is  the  wavenumber,  and  the  symbols  q  and  s  are  given  by 


5 


k2  -  k2  =  k2[l  -  (c/c,)2] 


s2  =  k-  -  k|  =  k2[l  -  (c/ Cg) 2 ]  . 


where  kd  and  ks  are  the  dilatational  and  shear  wave  numbers,  respectively 

kd  =  -[ps/(\  +2G)P  =  k(c/cd) 


kg  =  u;(ps/G)^  =  k(c/cg),  (16; 

where  \  is  the  first  Lamb  constant,  G  is  the  shear  modulus,  and  ps  is  the 
density  of  the  plate  material.  The  subscripts  a  and  s  in  the  amplitudes  As , 
Ba,  Ds,  and  Ca  refer  to  the  two  possible  types  of  Lamb  waves,  antisymmetric 
and  symmetric.  The  displacement  components  u  and  w  are  derived  from  the 
potentials  by 


3x  3z 


W  =  if  +  ^  *  (17) 

The  amplitudes  are  related  through  the  boundary  conditions,  which  are 
the  values  for  the  normal  and  shear  stresses  applied  at  the  two  faces  of  the 
plate.  These  stresses  are  found  in  terms  of  the  amplitudes  by  means  of  the 
following  equation 

Jx  =  \e  +  2Gcx 
cz  =  \e  +  2Gcz 

“zx  =  Gtzx  •  (18) 

The  elements  of  the  strain  tensor  are  given  in  terms  of  the  displacement 
vector  of  a  particle  with  components  u,  v,  w  by 


JW  _JU 

:zx  “  3x  +  3z 


This  leads  to  the  following  expressions  for  the  stresses,  using  Eqs.  (14) 
and  (17)  (the  factor  exp[i(.»t-kx)  ]  is  suppressed) 


(20) 


oz  =  G[As(k“  +  s~)cosh(qz)  +  Ba(k-  +  s2)  sinh(qz) 

-  Ca  2iks  sinh(sz)  -  Dg  2iks  cosh(sz)] 

Jzx  =  -G[As  2ikq  sinh(qz)  +  2ikq  cosh(qz) 

+  C  (k“  +s“)cosh(sz)  +  D,(k“  +s2)sinh(sz) ] .  (21) 

3  S 

The  applied  stresses  for  a  plate  in  vacuum  are  zero.  In  that  case  the  anti¬ 
symmetric  and  symmetric  waves  independently  satisfy  the  zero  boundary  condi¬ 
tions,  and  thus 

Ba(k-  +s-)  sinh(qd)  -  Ca  2iks  sinh(sd)  =  0 
Ba  2ikq  cosh(qd)  +  Ca(k2  +  s2)cosh(sd)  =  0 
As  (k“  +  s”)cosh(qd)  -  Dg  2iks  cosh(sd)  =  0 

and 

As  2ikq  sinh(qd)  +  Dg(k"  +  s“)sinh(sd)  =  0.  (22) 

By  elimination  of  the  amplitudes  from  these  equations,  one  arrives  at  the 
dispersion  relations 

(k2  +s2)2  sinh(qd)cosh(sd)  -  4k-qs  cosh (qd) sinh(j d)  =  0  (23) 

for  antisymmetric  waves,  and 

(k2  +s2)2Cosh(qd)sinh(sd)  -  4k2qs  sinh(qd)cosh(sd)  =  0  (24) 

for  symmetric  waves.  At  large  values  of  kd  both  dispersion  relationships 
approach  the  same  relation, 

4k2qs  -  (k2  +  s2)2  =  0,  (25) 

which  is  the  dispersion  relation  for  Rayleigh  waves  defined  as  waves  at 
the  surface  of  a  semi- inf inite  solid. 

The  space  dependence  of  the  (partial)  pressures  in  the  fluids  is  given 
by  the  following  expressions.  For  the  pressure  p^  of  the  incident  wave, 

Pi  =  exp[-ikG(x  sine  -  zcosc)  ]  .  (26) 

where  kQ  is  the  wave  number  of  the  wave  in  the  fluid.  Fox  the  pressure  pr  of 
the  reflected  wave, 

pr  =  Pr  exp[-ikQ(x  sine  +  zcoso)].  (27) 

For  the  pressure  p  of  the  transmitted  wave, 

Pt  =  Pt  exp[-ik'(x  sine'  -  zcosf')],  (23) 

where  k'  is  the  wave  number  of  the  transmitted  wave. 


The  angles  9  and  5'  are  che  angles  of  the  incident  and  transmitted  rays  with 
the  normal  to  the  face  of  the  plate,  respectively.  The  total  pressure  in 
the  fluid  at  the  surface  of  the  plate  on  the  side  of  the  incident  wave  is 

p(z  -  d)  =  pi  +  pr  =  (PQ  +  P)exp(-ikx),  (29) 

where  PQ  =  P^  exp(ikQd  cos:)  and  P  =  Pr  exp(-ik0d  cos:-).  The  pressure  at 
the  surface  of  the  plate  on  the  opposite  side  is 

p(z  =  -d)  =  P'exp(-ikx),  (30) 

where  P'  =  P  exp  (-ik'd  cos:)').  Here  the  coincidence  condition  is  used, 
namely 


k  =  kQsino  =  k'sint'.  (31) 

The  following  boundary  conditions  [Eqs.  (32)  through  (38)]  are  applicable 
to  the  problem  of  wave  propagation  in  a  fluid-loaded  plate.  At  the  side 
of  the  incoming  wave,  where  z  =  d,  one  has  continuity  of  the  normal  stress 
:z  and  the  shear  stress  cZK  =  0.  Thus,  combining  Eq .  (20)  with  Eq ,  (29), 
and  using  Eq.  (21),  one  has 

G [ As (k“+s-) cosh (qd)  +  Ba (k^+s2) sinh (qd) 

-  Ca  2iks  sinh (sd)  -  Ds  2iks  cosh(sd) ]  =  -  (PQ  +  P)  (32) 


and 


A„  2ikq  sinh(qd)  +  2ikq  cosh(qd) 

o  a 

+  C  (k^+s^)cosh(sd)  +  D  (k-+s^) sinh(sd)  =  0.  (33) 

a  S 

in  a  similar  fashion,  one  has  two  conditions  for  the  stress  at  the  opposite 
face,  where  z  =  -d,bv  combining  Eq.  (20)  with  Eq .  (30),  and  using  Eq.  (21) 

G[Ag (k“+s^)cosh(qd)  -  Ba (k2+s2) sinh (qd) 

+  Ca  2iks  sinh(sd)  -  Ds  2iks  cosh(sd)]  =  -  P'  (34) 

and 


-  As  2ikq  sinh(qd)  +  Ba  2ikq  cosh(qd) 

+  C  (k“+s^) cosh (sd)  -  D  (k2+s2)sinh(sd)  =  0.  (35) 

3  s 

Also,  at  the  interfaces,  the  component  of  the  particle  velocity  in  the 
solid  in  the  z  direction  has  to  equal  that  in  the  fluid.  This  leads  to  the 
relation  between  pressure  gradient  in  the  fluid  and  particle  acceleration 
in  the  solid  at  a  boundary.  In  general,  one  has 


o 


(36) 


The  set  of  boundary  conditions  can  be  advantageously  represented  in  the 
form  of  a  matrix  of  the  coefficients,  given  in  Table  I.  The  rows  are 
indicated  by  the  numbers  of  the  equations,  and  the  columns  are  marked  by 
the  corresponding  dimensionless  amplitudes  for  the  displacements  and  pres¬ 
sures.  The  density  applies  to  the  fluid  on  the  insonified  side,  - '  is 
the  density  of  the  fluid  on  the  other  side  of  the  plate.  The  relation 
c-  =  G/ps  has  been  used  to  arrive  at  the  given  form  of  the  matrix.  Since 
the  usual  FORTRAN  library  routine  for  complex  numbers  carries  only  circular 
functions,  a  transformation  is  made  whereby  s  =  is'  and  q  =  iq'.  Then 
the  hyperbolic  functions  are  transformed  into  circular  functions.  The 
result  is  shown  in  Table  II. 

To  find  the  dispersion  relation  for  free  waves  in  a  plate  loaded  by 
fluid  on  one  side  only,  one  omits  from  the  matrix  the  sixth  row  and  the 
fifth  and  seventh  columns  and  sets  the  determinant  value  of  the  remaining 
matrix  equal  to  zero.  The  dispersion  relation  for  waves  in  a  plate  loaded 
on  one  side  is  given  by 

iiais  +  (po/ps)q'd(ksd)4Uasin(q'd)sin(s'd)-Ascos(q'd)cos(s'd)] 

/ [ 2 (kd) cot6  ]  .  (39) 

If  one  indicates  the  matrix  elements  of  Table  II  by  m..,  the  symbols  ia,ls 
are  defined  by  ia  =  m^a^  -  m12m21’  anc*  ds  =  m33m44  “  mj4m43- 

The  dispersion  relation  for  a  plate  loaded  on  both  sides  by  the  same 
fluid  is  obtained  hy  omitting  column  five  from  the  matrix  and  setting  the 
determinant  value  of  the  remaining  matrix  equal  to  zero.  The  density  p' 
is  replaced  by  o0  and  d'  =9.  This  determinant  can  be  factored  in  two 
parts,  corresponding  to  antisymmetric  waves  and  symmetric  waves.  The  two 
dispersion  relations  are 

i-a4s  ~  (Po/^s)q'd(ksd)4Asc°s(q'd)cos(s'd)/[(kd)cote]  =  o  (40;) 

for  antisymmetric  waves,  and 
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Table  I.  Matrix  of  coefficients  of  equations  describing  wave  propagation 
in  fluid-loaded  plate  in  terms  of  hyperbolic  functions. 

Exact  elasticity  theory. 

Field  variables  - 


3  /J~ 

a 

0  /d  -  \  'd  - 

a  •» 

T)  'd 

■> 

P  ■'<:  o  ) 

.>  os 

?/  .  -  '  I 

O  S 

PV  (a1..;; 

( k*+s~  )d  sinh(qd) 

-2 iksd  “Slnh<  sd)  Kl 'cosh<  qd  > 

-J  iki^J  *  v- o s h <  sd  ) 

/ 

J  '  s 

0  *  s 

) 

2ikqd~cosh(qd) 

(k‘+-s~)d*cosh(sd>  2  ikqd  *  sinhi  qd) 

'*.  +%■  'd‘  sinh<  sd) 

) 

■) 

-( k"+s* >d -sinh(qd) 

2 iksd - sinhi sd )  ( k *  ”) d *co  sht qd ) 

-2  iksd  - COS.i(  ad  ) 

- 

) 

-  * ;  J  . 

s 

-2 ikqu-cosni qd) 

-<  k~  '»d  *  doshi  sd  )  2  ikqd  "s  mb'  qd ) 

•'  k*+s“  )d  -  sinht  sd) 

0 

0 

n 

’d  c  o  s  h  (  qd ) 

-ikd  coshvsd^  qd  sinhtqd) 

-ikd  sinh(sd) 

-ikd  cut; 

ikd  cot; 

<ksdi: 

(ksd  V 

qd  coshiqd) 

-ikd  coshlsd)  -qd  sinh»qd> 

|  ikd  s  tnh(  sd) 

) 

ikd  ' 

fk,d)- 

Table  II.  Matrix  of  coefficients  of  equations  describing  wave  propagation 
in  fluid-loaded  plate  in  terms  of  circular  functions. 

Exact  elasticity  theory. 

Field  variables-* 

UJiz  cji'-  A s/d2  lDa/d: 


— 

( k--s‘  ‘•)d^sin(q?d) 

tH,  ,  , 
2iks  d-sin(s  d) 

[k*-s*  t)d‘co8(q*d) 

. 

-2  Lks’d^cosCsM) 

eo/e. 

:  /  ; 

0  s 

-2ikq'd2cos(q'd) 

- (k‘ -s ' * )dfcCQ8<S 'c 

H2ikqfd‘ sin(q*d) 

-  -- 

-( k--s* 2)d29in(s’d) 

0 

0 

0 

-- 

-(k2-s* ‘)d2sin(q'd) 

-2iks 'd23in(s ’d) 

'k2-«*  2)d2cos(q'd) 

-2iks,d2cos(s,d) 

0 

. 

1 

■^Zlkq'd^cosiq’d) 

•*kz-s’J)<i2cos(s'd) 

•♦■ZikqM^sinCq'd) 

-  ( k  -  -  s 1 ’)d :sin(a’d) 

-  -  -  -  - 

. 

0 

r  ; 

< 

0 

q’d  cos(q'd) 

-ikd  co s (s *d) 

-q’d  ain(q*d) 

-ikd  sln(s‘d) 

mu 

n 

0 

q*d  co  s  (  q  *  d ) 

-ikd  coa(s'd) 

q'd  Sln(q'd) 

ikd  3in(s*d) 

°  ! 

i  ■ 

0 

ikd  cote' 

fk,j) ; 
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(41) 


i4a-ls  +  (pQ/ps)q 'd(ksd)^Aasin(q 'd)  sin(s'd)/ [  (kd)cot0  ]  =  0 

for  symmetric  waves.  In  a  strict  sense  it  is  not  possible  to  have  the 
two  wave  types  separately  in  a  plate  loaded  on  one  side  only,  since  the 
boundary  conditions  are  not  satisfied  in  that  case.  For  comparison,  one 
might  still  force  the  situation  by  altogether  leaving  out  rows  and  columns 
belonging  to  one  wave  type  or  Che  other.  That  way  two  more  dispersion 
relations  are  found  for  a  plate  loaded  by  a  fluid  on  one  side,  namely 

iAaas  -  (p0/ps)q'd(ksd)4Agcos(q'd)cos(s'd)/[2(kd)cos6)  =  0  (42) 

for  antisymmetric  waves  only,  and 

iAaAs  +  (p0/ps)q'd(ksd)4  Aasin (q ' d) sin (s ' d) / [2 (kd) cos0 ]  =  0  (43) 

for  symmetric  waves  only. 

Dispersion  Relations  From  Thick-Plate  Theory 

The  description  of  thick-plate  theory  for  straight-crested  waves  in 
fluid-loaded  plates  is  given  in  references  6  and  7.  The  matrix  of  the  coef¬ 
ficients  of  the  equations  of  motion  and  boundary  conditions  is  given  in 
Table  III.  The  symbol  y  indicates  the  phase  speed  c  of  the  waves  divided 
by  the  shear  wave  speed  cs.  The  subscripts  p  and  d  indicate  extensional 
and  dilatational  waves,  respectively.  The  correction  factor  for  the  ef¬ 
fective  shear  modulus  is  for  flexural  waves  and  K9  for  extensional  waves. 
U  is  the  average  displacement  in  the  z  direction  in  extensional  waves,  W 
is  the  average  displacement  in  the  x  direction  in  flexural  waves,  x  is  the 
average  strain  in  the  z  direction  in  extensional  waves,  and  <j>x  is  the 
average  angle  of  rotation  of  a  cross  section  in  flexural  waves. 

The  various  possibilities  of  fluid  loading  are  parallel  to  those  in 
exact  elasticity  theory.  The  dispersion  relation  for  a  plate  loaded  by 
a  fluid  on  one  side  is  obtained  from  the  matrix  in  Table  III  by  omitting 
columns  5  and  7  and  row  6.  The  determinant  value  of  the  resulting  matrix 
is  equated  to  zero  resulting  in 

i4aAs  "  (kgd)  (p0/ps)  (Aa  n33  +  n^/cotS  =  0.  (44) 

The  matrix  elements  are  indicated  by  n^j,  and  further  Aa  =  n^n01  -  n^n-,^, 
and  =  n33n44  -  n-^n^.  Again,  it  is  not  theoretically  possiSle  to“ 
satisfy  the  boundary  conditions  for  a  plate  loaded  by  a  fluid  on  one 
side  by  antisymmetric  wave  or  antisymmetric  waves  separately.  For  compari¬ 
son,  one  might  force  the  issue  by  just  leaving  out  the  pertinent  lines  from 
the  matrix.  This  leads  to  the  following  dispersion  relations.  For  anti¬ 
symmetric  waves,  fluid  on  one  side  only  one  has 

14a  -  !SY3(k3d)  (P0/pg)  n^/cotd  =  0,  (45) 

and  for  symmetric  waves,  fluid  on  one  side  only  has 

14g  -  •iY^(kgd)  (pQ/ps)  n33/cot6  *  0.  (46) 
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Table  III  -  Matrix  of  coefficients  of  equations  describing 
wave  propagation  by  fluid-loaded  plate  in 
thick-plate  theory 


Field  variables  -*• 

:x  w/d  l'/d  x  P;/(tDcs2)  P/<=0cs2)  P'/(;0cs:) 


-joeirc-,2- 

2 

1( kd)  <  j 2 

0 

0 

0 

0 

0 

i(kd)<,  2 

(kd)- 

0 

2 

c  llz 
o  s 

o  / 2p 

0  s 

-p  /2p 

0  s 

0 

0 

(kd)2(Yd2-Y2) 

. 

i<kdXvd2-2) 

0 

0 

0 

0 

° 

-i(kd)(yd:-2) 

m 

P  /2c 
o  s 

o  /2s 

o  s 

_ 

p0/:ps 

0 

v2(kd) 2 

0 

Y:(kd): 

-i(kd)cotd 

i(kd)cot? 

0 

— 

0 

Y  2(kd) 2 

0 

-Y:(kd)2 

0 

Q 

-i(kd)  ( o  /oMcot*' 
o 

If  the  wave  is  loaded  by  the  same  fluid  on  both  sides,  the  corresponding 
matrix  is  found  from  Table  in  by  omitting  column  5  and  setting  *  o'  and 
r  =  9 ' .  Then  the  determinant  value  of  the  matrix  can  be  split  in  two  fac¬ 
tors  for  antisymmetric  and  symmetric  waves  separately — with  the  dispersion 
relations,  for  antisymmetric  waves 

iA^  -  y3(kgd)  (po/ps)nn/cot9  =  0,  (47) 

and  for  symmetric  waves 

iAg  -  Y3Osd)  (p0/pg)n33/cote  *=  0.  (48) 


Discussion  of  Results 

In  this  section  results  are  discussed  of  the  application  of  the  root 
finder  to  wave  propagation  in  plates.  Only  those  complex  roots  are  given 
that  originate  in  zero  order  antisymmetric  and  symmetric  Lamb  waves  in  an 


unloaded  plate.  The  case  of  steel  was  chosen,  with  Poisson's  ratio  equal 
to  0.3028,  a  shear  wave  speed  cs  =  3264  m/s,  and  water  of  4°C  temperature 
for  which  c0  =  1447  m/s.  The  value  of  the  Rayleigh  wave  speed  relative 
to  the  shear  wave  speed  in  steel  is  0.9278,  and  this  was  chosen  as  the 
value  of  the  correction  factor  for  the  effective  shear  modulus  in  the  thick- 
plate  computations.  The  results  for  the  various  options  are  shown  in 
Tables  IV  through  VII.  The  results  of  the  change  in  the  real  part  of  the 
relative  wave  speed  Rey  are  not  overly  significant.  First,  the  variation 
with  frequency  is  partly  due  to  the  lack  of  accuracy  of  the  real  seed  used 
in  the  program;  thus  no  trends  are  visible.  Secondly,  the  fractional  change 
in  Rey  is  very  small.  More  important  is  the  behavior  of  the  imaginary  part 
of  y.  Two  aspects  of  the  results  will  be  discussed:  in  the  first  place, 
just  how  noticeable  is  the  effect  of  the  presence  of  two  wave  types,  sym¬ 
metric  and  antisymmetric;  and  in  the  second  place,  how  close  are  the  ap¬ 
proximations  of  the  thick-plate  theory  to  the  results  of  exact  elasticity 
theory? 

In  Figs.  2  and  3  a  comparison  is  made  between  the  results  from  exact 
elasticity  theory  for  a  plate  loaded  on  one  side  by  a  fluid.  In  Fig.  2  the 
complex  root  originates  in  an  antisymmetric  wave  in  an  unloaded  plate.  In 
Fig.  3  it  derives  from  a  symmetric  wave  in  an  unloaded  plate.  In  both 
cases.  Curve  if  1  represents  the  result  of  admitting  both  antisymmetric  and 
symmetric  waves  to  propagate  in  the  plate.  This  is  the  correct  way  of 
satisfying  the  boundary  conditions  on  each  side  of  the  plate.  If  one  leaves 
out  the  second  type  of  wave,  the  symmetric  one  in  Fig.  1  and  the  antisym¬ 
metric  one  in  Fig.  2,  Curves  if 2  are  obtained.  The  curves  show  that  admit¬ 
ting  both  types  of  waves  instead  of  only  one  gives  a  very  small  difference 
at  lower  values  of  the  dimensionless  wave  number  kd,  but  the  difference 
becomes  quite  pronounced  at  higher  values  of  kd,  above  kd  s  6.  At  this 
value  of  kd,  the  wavelength  is  of  the  order  of  the  thickness  of  the  plate. 

The  comparison  of  exact  theory  with  thick-plate  theory  in  Figs.  4  and  5 
shows  a  corresponding  behavior: the  thick-plate  theory  gives  results  reason¬ 
ably  close  to  exact-plate  theory,  except  in  the  region  where  the  exact  theory 
starts  to  display  the  influence  of  the  complementary  wave  type,  around 
kd  s  6  or  7.  Thus,  even  where  the  thick-plate  theory  admits  both  types 
of  waves  to  take  part  in  the  propagation,  it  is  not  able  to  correctly  pre¬ 
dict  the  strong  reduction  in  attenuation  that  occurs  at  higher  values  of 
kd  according  to  the  more  reliable  exact  theory.  In  thick-plate  theory  there 
is  very  little  change  in  y  when  both  types  of  waves  are  admitted  as  compared 
with  one  type  only. 
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Table  IV.  Roots  of  dispersion  relation  in  thick-plate  theory  corresponding  to 
antisymmetric  waves  in  unloaded  plate. 


UNLOADED  PLATE 

LOADED  PLATE 

A+S.  loaded  one  side 

A  only,  loaded  one  side 

A  onlv,  loaded  both  sides 

> 

Re (Ay) 

Im(Ay) 

Re (Ay) 

lm(  '.y) 

ReC *>) 

Im  (  '  ! ) 

0.7 

0.5322 

0.37254 

3.04x10-* 

♦3x10-6 

2.0767x10-2 

1 3  x 1 0~  6 

1.69x10-* 

t3xl0-6 

2.0775«10-2 

t3»10-6 

1.72x10"* 

t7xl0"6 

4.04  52«10-2 
t7xl0-b 

1.1 

0.67  58 

0.74338 

1. 0965x10-2 
t3xl0“6 

1.83x10-* 

t3xl0~6 

1.0964xl0-2 

t3»10-6 

-7 ,25*10"* 
i5xl0“6 

2.1885-10-2 

t5xlO-f 

1.5 

0.7575 

1.13625 

1.9«10-5 

two-6 

8 . 24  3  x  L0" 3 

two-6 

-9.1 xl0‘5 
tlxio-6 

8.233  «10”J 
tlxio-6 

-3.79x10-* 

t3«10"6 

1.64  5  6-10-2 
t3xlO"6 

2.5 

0.84  96 

2.12400 

1.18x10-* 

tWO-6 

5.472x10'’ 

two-6 

1 .2  xlO"5 

two-6 

5.473  xlO"3 
tl«10-b 

1 . 0946-10-2 
±3xl0-6 

4.5 

0. 90Uj 

4.05135 

1 . 568x10"* 
t6xl0"7 

3. 3481* 10“ * 

i6*10“7 

9.00X10-6 

t’xlO-7 

3 ,3512*10” * 

3. 5  xlO-5 
tl-10"6 

6.699x  10"1 
*3xi  -j-o 

b.  5 

0.9142 

5.94230 

1.07x10“* 

two-6 

2.392  xlO"’ 

tlxio-6 

-4.72X10-5 

tlxlO-7 

2. 3979x10-* 
tlxlO-’ 

-6. 9x10“ 5 
tl-10-6 

4.795x10-’ 

tW0“6 

8.5 

0.9197 

7.31745 

1.611x10-* 

t5»10"7 

1.8515x10"- 

t5«10“7 

5.1*10“6 
t  2  * 1 0“ 7 

1.86*10"3 

t2y10“~ 

-8.i  no-6 

t2*10"' 

3.722-10-3 

t2«10-' 

I 

10. 1 

0.9220 

1.393x10-* 

t5xl0'7 

1. 5651*10“  ’ 
i5*10“7 

:.:4*io-5 

Tft*10”7 

3 .1 501 • 10"  j 
*6*10“ " 

Steel,  v  *  0.3028  Water,  c0  ■  1447  m/9 

c3  •  3264  m/s  y  R  In  steel  »  0.9278 


14 


,1' 


Table  V.  Roots  of  dispersion  relation  in  exact  elasticity  theory  corresponding 
to  antisymmetric  waves  in  unloaded  plate. 


UNLOADED  PLATE 


0.5338 

0.37366 

0.6801 

0.74811 

0.7643 

1.14645 

0.8601 

2.15050 

0.9122 

4.104  9  0 

0.9240 

6.00600 

0.9269 

7.87865 

0.9267 

9.368  7  6 

LOADED  PLATE 


A+S,  loaded  one  side  A  only,  loaded  one  side!  A  onlv,  loaded  NHh  side 


Re(Ay)  Iih(Ay)  Re(^y)  ImCA-y)  Re(<lv) 


tl*10-6  ±lxlO-6 


8.88xl0"5 

±5xL0-7 


±5xl0-7  ±5«10-7 


±4<10"7  ±4xlQ-7  I  ±1«10-7 


4.0471x10-3  |  1 ,198x 10““  4.1106xio-3 

±2“10-7  :2xl0-7 


2.56»10*S  8, 7082x10* 3 

±2xl0"5  ±9x10" 7 


±2xl0"6  ±2xl0"6  ±2«10"5  ±2«10"6  ±2*10"5 


a.SSxlO"1’  9.57x10'*  4.549x10" 

4x10-3  rlxlO"6  tlxlO-6 


8.2174«10-3 

±5x10-' 


Steel,  v  •  0,3028  Water,  c0  ■  1447  m/ s 

c  »  3264  m/ s  yD  In  steel  -  0.9278 


Table  VII.  Complex  roots  of  dispersion  relation  in  exact  elasticity  theory 
corresponding  to  symmetric  waves  in  unloaded  plate. 


Im  (y) 


Curve  #1  -  Antisymmetric  and  symmetric  wave 
occurring  simultaneously. 

Curve  If 2  -  Antisymmetric  wave  only. 


Fig.  2.  Imaginary  part  of  relative  phase  speed  y  as  a  function 
of  the  dimensionless  wave  number  kd,  corresponding  to 
antisymmetric  root  of  the  unloaded  case.  Exact 
elasticity  theory. 
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Curve  #1  -  Symmetric  and  antisymmetric  waves 
occurring  simultaneously. 

Curve  #2  -  Symmetric  wave  only. 

Imaginary  part  of  relative  wave  speed  y  as  a  function 
of  dimensionless  wave  number  kd,  corresponding  to 
symmetric  root  of  the  unloaded  case. 

Exact  elasticity  theory. 


Fig.  4.  Comparison  of  imaginary  part  of  relative  phase  speed  y 
as  a  function  of  dimensionless  wave  number  kd  for 
exact  elasticity  theory  (Curve  // 1 )  and  thick-plate 
theory  (Curve  #2) ,  for  root  corresponding  to 
antisymmetric  root  in  unloaded  case. 
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APPENDIX  A 


DESCRIPTION  OF  ROOT-FINDING  PROGRAM 


The  principle  of  the  complex  root  finder  described  in  this  report  is 
general,  but  the  specific  form  of  the  program  is  adapted  to  the  problem  of 
finding  roots  of  the  dispersion  relation  of  straight-crested  waves  in  a 
fluid-loaded  plate.  The  development  of  the  dispersion  relation  for  such 
waves  is  given  in  the  Application  section  of  this  report.  Flow  diagram, 
list  of  symbols,  source  listing,  and  examples  are  attached  as  Appendices  B, 

C,  D,  and  E,  respectively. 

First  the  material  parameters  are  entered  into  the  program:  POI, 
Poisson's  ratio  of  the  plate  material;  COCS,  the  ratio  of  the  sound  speed 
in  the  fluid  to  the  shear  wave  speed  in  the  plate;  and  RH>1,  the  ratio  of 
the  fluid  density  to  the  density  of  the  plate  material.  Next,  the  complex 
correction  factor  (AKAR,  AKAI)  for  the  effective  shear  modulus  in  thick- 
plate  theory  is  entered,  with  the  (uncritical)  value  of  the  step  in  z,  DELZ, 
used  to  approximate  the  derivative  3f/3z  (see  Complex  Root-Finder  Program 
section).  Then  the  frequency  is  entered  in  the  form  of  the  dimensionless 
wave  number  kgd  (AKSD)  for  shear  waves  in  the  plate,  with  the  corresponding 
dimensionless  propagation  speed  for  an  unloaded  plate,  ZNOT,  which  should  be 
computed  by  a  program  for  finding  real  roots  [8].  On  the  same  line,  a  maxi¬ 
mum  number  of  iterations  are  entered,  which  serves  to  provide  an  exit  from 
the  program  in  case  it  fails  to  converge.  Three  control-type  variables  are 
entered  next:  ANR,  the  number  of  steps  in  which  the  approach  from  zero 
density  to  full  density  is  performed;  TOL,  the  number  for  the  relative  ac¬ 
curacy  of  the  final  result;  and  the  multiplication  factor  FSTEP ,  for  the 
step  size  in  the  movement  of  the  test  pairs  (see  Complex  Root-Finder  Pro¬ 
gram  section).  The  option  number  BID  chooses  from  various  options,  concern¬ 
ing  exact  theory  or  thick-plate  theory  and  antisymmetric  cr  symmetric  waves 
and  fluid-loading  on  one  side  or  both  sides.  The  listing  of  these  options 
is  given  in  Appendix  B. 

The  computations  and  loops  follow  closely  the  discussion  of  the  main 
algorithm  in  the  Complex  Root-Finder  Program  section.  If  the  program  fails 
to  converge,  a  statement  is  printed  out  which  indicates  an  exit  left  or 
right  or  up  or  down,  depending  on  the  direction  of  motion  of  the  pertinent 
test  pair.  The  name  SIGN  refers  to  the  quantity,  the  algebraic  sign  of  which 
is  the  criterion  for  the  location  of  the  root  with  respect  to  the  test  pair. 
If  SIGN  is  equal  to  zero,  a  message  to  that  extent  is  printed  and  the  program 
transfers  to  the  end  where  a  return  option  number  is  to  be  entered. 

The  accuracy  of  the  final  answer  is  determined  by  the  comparison  of  the 
larger  of  the  relative  errors  in  the  real  and  imaginary  parts  of  the  change 
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FuECEOlNG  PAGE  BLANK -WOT  Flu»iD 


in  reiacive  wave  speed  co  the  number  TOL  entered  in  the  program.  For  every 
step  in  the  density,  an  intermediate  value  of  the  complex  wave  speed  is 
printed  out.  A  dimensionless  attenuation  factor  -(kid;  is  computed  and 
printed  out  that  is  computed  according  tc 


(ksd)vi 

-Odd)  =  — - v  . 

YI  +  y2 


( A 1  y 


where  1,2  refer  to  the  real  and  imaginary  parts  of  a  complex  number,  ks  is 
the  wave  number  for  shear  waves,  2d  the  thickness  of  che  place,  yj  and  Y2 
are  che  real  and  imaginary  parts  of  the  dimensionless  propagation  speed 
'<  -  c/cs,  c  is  the  propagation  speed,  and  cs  is  the  speed  of  shear  waves  in 
the  place.  Ac  the  end  of  the  program,  a  return  option  number  is  requested 
chac  returns  the  program  to  any  of  the  five  lines  where  input  data  are 
entered.  Typing  in  zero  for  the  option  number  results  in  an  exit  from  the 
program. 


Special  caution  is  due  in  the  use  of  a  complex  square-root  routine, 
since  the  square  root  is  a  double-vaiued  function.  It  occurs  several  times 
in  che  program:  first  in  computation  of  the  quantities  q'  and  s',  defined  by 
q'-  =  kg  -  k- ,  and  s'2  =  k|  -  k2  [compare  Eqs.  (14)  and  (15)].  The 
two  possible  roots  differ  by  a  factor  of  -1,  of  course.  The  structure  of 
the  dispersion  determinant  is  such  that  che  change  from  one  root  to  the  ocher 
introduces  a  factor  of  -1  into  che  determinant  value,  and  thus  does  not  affec 
che  location  of  zeros  of  the  determinant.  The  program  converges  properly  as 
long  as  one  stays  with  one  branch  of  the  function  and  does  not  jump  from  one 
branch  to  the  other.  The  FORTRAN  subroutine  selects  always  the  square  root 
that  has  a  positive  real  part;  if  the  real  part  is  zero,  it  selects  the 
positive  purely  imaginary  root.  This  amounts  co  a  branch  cut  along  the 
negative  real  axis  in  the  complex  plane  of  the  argument  of  the  function,  and 
this  prevented  convergence  of  the  program  in  some  cases.  The  problem  was 
remedied  by  the  following  method.  If  the  numerical  value  for  q ' 2  or  s'“  has 
a  negative  real  part,  it  is  multiplied  by  -1  before  the  square  root  is  taken. 
The  result  is  then  multiplied  by  -i.  This  is  equivalent  to  moving  the 
branch  cut  from  the  negative  real  axis  to  the  positive  real  axis.  This 
should  solve  the  problem  in  all  cases  where  the  root  is  not  close  to  the 
origin,  as  compared  with  the  chosen  stepsize. 


In  cne  second  place  the  square-root  routine  appears  in  the  calculation 
of  cot:,  where  5  is  Che  angle  with  Che  normal  to  the  plate  of  the  radiation  . 
emitted  into  the  fluid.  It  is  calculated  by  che  equation  c ocf  =  ( (c/c0) “-1 j  : 
which  follows  from  che  coincidence  condition  between  the  waves  in  the  plate 
and  in  the  fluid.  Since  the  imaginary  part  of  the  phase  speed  c  is  always 
small  compared  with  the  real  part,  there  is  no  problem  in  this  case  with  the 
branch  cut  of  the  subroutine,  as  long  as  c  is  larger  than  cQ  and  not  very 
close  to  c0.  If  c<c0,  one  has  che  situation  where  there  is  no  raaiucion  into 
the  fluid;  the  factor  [  (c/c.3j  —  1  ] *-  is  purely  imaginary,  and  the  root  finder 
is  not  aole  to  operate.  Therefore  the  program  checks  for  the  aigeDraic  sign 
of  the  expression  (c/c0)~-L  in  the  beginning,  wnen  c  has  still  the  value  of 
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che  phase  speed  for  an  unloaded  plate.  If  this  is  negative,  the  program 
prints  a  statement  chat  no  radiation  into  the  fluid  is  possible  and  moves 
to  che  request  for  a  new  return  option  number. 
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COMPUTE 


28 


'REL . 
ERROR 
-  TOL 


FV1  *  FV2 

V' 


CALCULATE 
REGA,  AMGA 


PRINT 

GAMMA 

-ZNOT, 

ERROR 


'  PRINT 
DIM.  LESS 
ATTEN.  i 
FACTOR  / 


30 


FLOW  DIAGRAM  FUNCTION  SUBPROGRAM  DSP 


33 


► 


FLOW  DIAGRAM  FUNCTION  SUBPROGRAM  DSP2 

(CALLED 

BY  OS  J 

zr 

SET  FF 
DEPENDING 
ON  BID 


I 


COMPUTE  DSP2 
FOR  GIVEN  BID 


APPENDIX  C 


AI 

AKAP 
AKAR 
AKAI 
AKFD 
AKSD 
ANR  = 
AMGA 
AMI  = 
ATF 
BID  = 


LISTING  OF  MAJOR  SYMBOLS  IN  FORTRAN  PROGRAM 

i,  imaginary  unit 

K,  complex  correction  factor  for  shear  modulus 

Real  part  of  k 

Imaginary  part  of  K 

k  d,  dimensionless  wave  number  in  fluid 
o 

k  d,  dimensionless  wave  number  for  shear  waves  in  plate 
s 

NR  Number  of  steps  to  reach  nominal  density  RHM 

Imaginary  part  of  wave  speed  y 

RAM  Number  of  iterations  allowed  in  loop 

-k^d,  where  k0  =  Im(k) ,  dimensionless  attenuation  factor 

ID  Number  of  required  option 

1:  Thick-plate  theory,  antisymmetric  and  symmetric  waves, 

fluid  on  one  side. 

2:  Thick-plate  theory,  antisymmetric  wave  only,  fluid 

on  one  side. 

3:  Thick-plate  theory,  antisymmetric  wave  only,  same 

fluid  on  both  sides. 

4:  Thick-plate  theory,  symmetric  wave  only,  fluid  on 

one  side. 

5:  Thick-plate  theory,  symmetric  wave  only,  same  fluid 

on  both  sides. 

6:  Exact  theory,  antisymmetric  and  symmetric  waves, 

fluid  on  one  side. 

7:  Exact  theory,  antisymmetric  wave  only,  fluid  on  one  side. 

8:  Exact  theory,  antisymmetric  wave  only,  same  fluid 

on  both  sides. 

9:  Exact  theory,  symmetric  wave  only,  fluid  on  one  side. 

10:  Exact  theory,  symmetric  wave  only,  same  fluid  on 

both  sides. 
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cocs 

DS  (Z) 

DSP (Z) 

DSP2 (Z) 

DFDZ 

DEL 

DELZ 

DELP 

DELH,  DELV 

DZ 

F 

FF 

F1.F2, 

F3,F4 

FVl ,FV2 , 
FH1.FH2 

FSTEP 

GAMD 

GAMP 

ID 

IDSP ,  IBD 
K 


KAM 


c0/cs»  propagation  speed  in  fluid  relative  to  shear  wave 
speed  in  plate. 

Function  subprogram  calling  DSP(Z)  or  DSP2(Z). 

Function  subprogram  for  dispersion  relation  in  thick-plate  theory. 
Function  subprogram  for  dispersion  relation  in  exact  theory. 
Approximation  for  3f/3z. 

Real  step  size. 

Az,  used  in  computation  of  approximation  to  3f/3z. 

Intermediate  symbol  for  step  size  DEL. 

Step  sizes  for  movement  of  test  pairs,  horizontally  and  vertically. 
Complex  step  size 
[  (c/cq)2  -  1 ]  2 

Factor  in  dispersion  relation  FF  =  1 ,  fluid  on  both  sides 

FF  =  2  ,  fluid  on  one  side 


f(Zl),  f(z2),  f(z3),  f(z4) 


Flags  to  indicate  direction  of  motion  of  test  pairs. 


Factor  to  adjust  step  size  DEL. 

y^,  phase  speed  of  dilitational  waves  relative  to  shear  wave  speed. 


Y  »  phase  speed  of  extensional  waves  relative  to  shear  wave  speed. 


BID 


Intermediate  option  numbers. 

Return  option  number 

=  1,  program  returns  to  entering  POI,  COCS,  RHM 

*  2,  program  returns  to  entering  AKAR,  AKAI,  DELZ 

3  3,  program  returns  to  entering  AKSD,  ZNOT,  AMI 

=  4,  program  returns  to  entering  ANR,  TOL,  FSTEP 

=  5,  program  returns  to  entering  BID 

*  0,  program  exits. 


-  AMI 
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KR,  KL, 

KUP.KDOWN  Counters  in  loops  to  check,  number  of  iterations. 

NR  =  ANR 

NRV  Counter  in  advancing  the  parameter  RHO. 

P01  Poisson's  ratio. 

REGA  Real  part  of  relative  wave  speed  y . 

RHM  Nominal  density  of  fluid  loading  the  plate,  divided  by 

density  of  plate  material. 

RHO  Stepwise  varied  value  of  relative  density,  varying  from 

zero  to  RHM. 

QPR  q'd  =  [ (k,d)“  -  (kd)/]'s,  dimensionless  wave  number. 

SIGN  Expression,  the  algebraic  sign  of  which  determines  location 

of  root  relative  to  test  pair. 

SPR  s'd  =  [(kgd)“  -  (kd)z]’2,  dimensionless  wave  number. 

TOL  Limit  for  relative  error  in  the  root. 

Z  Variable  for  root  used  in  calling  function  subroutines. 

Z1,Z2  Vertical  test  pair. 

Z3,Z4  Horizontal  test  pair. 

ZNOT  Value  for  the  real  root  of  the  dispersion  relation  without 

fluid  loading,  from  a  real  root  finder,  serving  as  the 
seed  for  starting  the  program. 

ZVAR  Intermediate  value  of  root,  serving  as  seed  for  the  next 

step  in  the  parameter  RHO. 
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APPENDIX  D 


SOURCE  PROGRAM  LISTING 


COMPLEX  Z,DSP,AI ,DFDZ,DZ,ZVAR,Z1 ,Z2,Z3,Z4*F1 
COMPLEX  A  K  A  P  ,  D  S  »  C  M  T'  L  X 

COMMON  RHO  >  AKSD , GAMP . G AMD  »  AK AR , AKA  I  ,  AKFD,  I  [i , 

I  WRITE  <5.30) 

30  FORMAT <  '  $  E  N  T  E  R  POI » COCSf RHM J  ') 

READ  (5.3?)  POI,CQCS,RHM 

39  FORMAT <3F 15,0) 

3  WRITE  <5» 31) 

31  FORMAT ( ' IENTER  AKAR»  AKA I ,  DELZ  !  '< 

READ  (5,39)  AKAR  ,  A<AI ,  DEi,Z 
AKAP®CMPL>  <  AKAR , AKA  I ) 

II  WRITE  (5,505 

50  FORMAT (' $ENTER  A K 3 D , Z N 0 T , J  OF  IT.!': 

READ  (5,39;  AKSD , ZNO T , AM  I 
IF  (ZNOT/COCS-1 . )  61,51,62 

61  WRITE  (5,63) 

63  FORMAT  ('  NO  RADIATION  INTO  FLUID'', 

GO  TO  100 

62  K  A  M  =  A  M I 

10  WRITE (5,45) 

•45  FORMAT  (  '  $  E  N  T  E  R  ANR,  TOL,  FSTEPJ  ' 

47  READ < 5 , 39 ) ANR , TOL , FSTEP 

12  WRITE  (5,33) 

33  FORMAT (  ' $ENTER  BID:  '  ) 

READ  (5,8)  BID 
3  FORMAT  ("15. 0) 

I D  =  B I D 

AI=CMF'LX  (  0  .  ,  1  .  ) 

GAMP  =  SQRT ( 2-/(1. -PO  I  )  ) 

G AMD- SORT  (2  .  :K  ( 1  ,-POI )/(  1  .  -2  .  *POI  ■  ) 

AKFD®  AKSD /’CUCS 
R  H  0  =  0  . 

NR=INT( ANR ) 

NR  V  =  Q 

C  COMPUTE  STEPSIZE  DEL 

ZVAR=ZNOT 
Z=ZVAR+DELZ 
D  F  D  Z  ®  D  3 ( Z ) /DELZ 
RHO=RHO+RHM/ANR 
NRV=NRU*3 
D Z  =  - D S ( ZVAR  '  / DFDZ 
X=ABS ( REAL ( OZ I ) 

Y  =  ABS  (  A  I  MAG  (  DZ 1  ) 

DEL® AMAX 1  <  >  Y  ;<F  STEP 

WRITE(5,3C1 )  DEL 

P«EC£DIi\((j  pai ar 


»F2,F3»F4,D5F2 
I  DSP 


BLWK-M T  F3ii.fi) 


301 


FORMAT ('  DEL  -  • £12.5 > 

DELF=DEL 

C  FIND  INITIAL  DIRECTION  l.  F  MOVEMENT  OF 

C  VERTICAL  TEST  FAIR 

70  DEL = DEL P 

64  Zl=ZVAR-AI#DEl 
Z2=ZVAR+AI*DEL 

65  F 1  =  D  S ( Z  1  ) 

F 2  =  D S ( 22  ) 

NR-0 

KL  =  0 
K  U  P  =  0 
K  D  0  U  N  =  0 

SIGN=AIMAG(F1  )*REAL(F2)-REAL(F1  )  *A  I  h  AG  •  2. 

I F ( S I GN ) 7 1 , 72  »  7 3 

71  FH 1  =  -  1  . 

GO  TO  74 

72  UR  I TE ( 5 ,  302  > 

302  FORMAT  ('  LEFT-RIGHT  SIGN  =  ZERO') 

GO  TO  100 

73  FH  1  =  1  . 

74  DELH=FH1*DEL 

C  LOOP  TO  FIND  VERTICAL  LINES  BRACKET ] NG  ROOT 

75  Z1=Z1+DELH 
Z2=Z2+DELH 
F 1 =DS ( Z 1 > 

F  2  =  D  S  (  Z  2 ) 

S I GN  =  A I MAG ( F 1 ) *REAL ( F2 ) -A  I  MAG < F2 ) #RE AL < F  1  ) 
IF  (SIGN;  81,72,83 
81  FH2=-1. 

KL=KL+ 1 

IF  (KAM-KL)  99,99,84 
99  WRITE (5, 3 10) 

310  FORMAT  (  '  EXIT  LEFT  '  ) 

GO  TO  100 

S3  FH2  =  1. 

KR  =  KR+ 1 

IF ( KAM-KR )  98 , 98  >84 
98  WRITE  (5,311) 

311  FORMAT (  '  EXIT  RIGHT  '  ) 

GO  TO  100 

84  IF  (FH1*FH2)  92,75,75 

c  find  initial  direction  of  movement 

C  OF  HORIZONTAL  TEST  PAIR 

92  Z4=(Zl+Z2)/2.+DEL 

KR  =  0 
KL  =  0 

Z3=(Z1+Z2)/ 2,-DEL 
F  3  =  D  S  (  Z  3  ) 

F  4  =  D  S  (  Z  4  ) 

3  I G  N  =  A I M  A  G  (  F  3  >  *  R£  A  L  (  F  4  )  -  A  I  M  A  G  (  F  4  i  *  R  E  A  L  <  F  3  ) 
IF  (SIGN)  101,102,103 


101 


F  V 1  =  1  . 

GO  TO  104 

102  WRITE  (5,304) 

304  FORMAT  (  '  UP-DOWN  SIGN  =  2 ERO 
GO  TO  100 

103  F  V 1  =  -  1  . 

104  DELV=FV1*0EL 

C  LOOP  TO  FIND  HORIZONTAL  LINES  BRa  rEIIND  ROHI 

105  Z4=Z4+AI4DELU 
Z3=Z3+AI*DEL  0 
F3=DS(Z3) 

F4*DS<  Z4  ) 

SIGN=AIMAG<F3)*REAL <F4 )-AIMAG(F4  >*REAL  < F3 ) 

IF  (SIGN)  111,102,113 
1 1 1  FU2  =1  . 

NUP  =  KUP+ 1 

IF(KAM-KUP)  97,97.114 
97  WRITE(5»305) 

305  FORMAT  ('  EXIT  UP') 

GO  TO  100 

113  FV2  =  - 1  . 

KD0WN=KD0WN+1 

IF  (KAh-KDOWN)  96,06,114 
96  WRI TE ( 5 »  306  ) 

306  FORMAT ( '  EXIT  DOWN' ) 

GO  TO  100 

114  IF  <FU1*FV2)  122,105,105 

122  REGA  =  REAL(Z1  ) 

AMG A  =  A I  MAG  <  Z3 ) 

K  U  P  -  0 
K  D  0  W  N  =  0 

133  AMGA= ABS ( AMGA ) 

DREG=ABS(REGA-ZNOT) 

C  CHECK  IF  REQUIRED  PRECISION  IS  REACHED 

IF  (DEL/AMINl  ( DREG , AMGA ) -TOL >  132,132,  131 

131  CONTINUE 
BEL= DEL/10. 

Zl»<Z3+Z4)/2.-AI*BEL 
Z2=(Z3+Z4)/2. +AI*DEL 
GO  TO  65 

132  REGA=RE AL  (Zl)  -DELH/2. 

AMGA  =  A  I  MAG ( Z3 ) -BELV/2. 

DELH  =  ABS( DELH/2.  ) 

C  PRINT  PRELIMINARY  RESULTS  FOR  INTERMEDIATE  RHU 

WRITE (5, 140)  RHO»REGA»AMGA»DELH 
140  FORMAT ( '  RHO , REGA , AMGA , DELH= ' , 2F 1 0 , 6 . 2E 1 5 . 5  ) 

C  CHECK  IF  FINAL  RHO  HAS  BEEN  REACHED 

IF  <  NR-NRY )  145, 145, 146 
146  ZVAR=REGA+AMGA*AI 

RHO=RHO+RHM/ANR 
NRV=NRV+1 
GO  TO  70 
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145  GO  TO  (701 ,702,802,703, 303. 704, 705, 805,?<V'  . 8  >6  •  fi 

701  WR I TE ( 5 , 7  1  1  ) 

GO  TO  707 

702  WRITE(5»712) 

GO  TO  707 


302 

WRITE ( 

5,812) 

GO  TO 

707 

703 

WRI  TE 

(5,713) 

GO  ro 

707 

803 

WRITE 

(5,813) 

GO  TO 

707 

704 

WRITE 

(5,714) 

GO  TO 

70  7 

705 

WR I TE ( 

5,715) 

GO  TO 

707 

805 

WRITE ( 

5,815) 

GO  TO 

70  7 

706 

WRITE 

(5,716) 

GO  TO 

707 

806 

WRITE 

\  5 , 8 1  A  ) 

707 

CONTINUE 

C  COMPUTE  ATTENUATION  FACTOR 

ATF=AKSD#AMUA/  (  REGA  KrO  f  AMGA*»  ;:  ; 

711  FORMAT  < '  THICK  PLATE,  A.  AND  S> .  ,  FLUID  ONE  SIDE') 

712  FORMAT (  '  THICK  PLATE,  A. ONLY,  FLUID  ONE  SIDE') 

312  FORMAT ( '  THICK  PLATE,  A. ONLY,  FLUID  DOTH  SIDE' 

713  FORMAT ( '  THICK  PLATE,  S.ONLY,  FLUID  ONE  EIDK 

313  FORMAT  < '  THICK  PLATE,  S.ONLY,  FLUID  RUTH  SIDES') 

714  FORMAT ( '  EXACT,  A.  AND  S.,  FLUID  ONE  SIDE') 

715  FORMAT ( '  EXACT,  A. ONLY,  FLUID  ONE  SIDE') 

815  FORMAT ( '  EXACT,  A  ONLY,  FLUID  BOTH  SIDES') 

716  FORMAT ('  EXACT,  S.  ONLY,  FLUID  ONE  SIDE') 

816  FORMAT ( '  EXACT,  S.  ONLY,  FLUID  BOTH  SIDES') 

171  FORMAT  (2F10.6) 

WRITE  (5,172)  AKSD , 2N0T 

172  FORMAT (  '  UALUES  OF  AKSD » ZNQT ,  ARE  ',2F10.6) 

WRITE  (5,173* 

173  FORMAT  ('  REAL  AND  IM.  PARTS  OF  GAMMA--ZNOT  AND  ERROR  A 
REGA=REGA-ZNOT 

WRITE  (5,174)REGA,AMGA,DELH 

174  FORMAT ( 3E 1 5 . 5 ) 

WR I TE (5,175)  ATF 

175  FORMAT  ('  DIMENSIONLESS  ATTEN. FACTOR  =  ',E12.5: 

100  WRITE  (5,159) 

159  FORMAT  ( ' 4ENTER  RETURN  OPTION  NUMBER  ') 

READ  (5,151)  K 
151  FORMAT  (14) 

GO  TO  (1,2,11,10,12)  K 
END 
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FUNCTION  D  S  <  Z  ) 

COMF'LEX  DS,DSP»Z»DSP2,AKAF 

COMMON  RHO , AKSD,GAMP, GAMD , AK AR  , AK A  I  *  ANFD  , 1 0 , I  DSP 
BID* ID 

IF (BID -5. 5)50*50*60 
DS=DSP ( Z ) 

GO  TO  70 
DS  =  DSP2<  Z> 

RETURN 

END 


10 


20 


30 

40 


FUNCTION  DSP(Z) 

COMPLEX  DSP , A I , AKAP 

COMPLEX  A1 1 ,  A22,Si  1  »S22,DELA,DELS,Z,F,CMF'LX,CSGRT 
COMMON  RHO, AKSD, GAMP, GAMD, AKAR.AKAI »  A  K  F  0  » ID, IDSP 
GO  TO  <7,6,*,  ,6,  i.)  ID 
FF= 1 .  '  ' 

GO  TO  7 
FF- 2  . 

CONTINUE 

AI=CMPLX<0. 0,1.0) 

AKAP=CMPLX( AKAR.AKAI ) 

A1 1=AKSD**2*<  Z**2-GAMP*#2>  -'3  .  -AKAP**2#Z**2 
A22=ANSD**2*< AKAP**2-Z**2 ) 

DELA  =  A1 1*A22  +  AKSD**2*Z#*2*AKAF'**4 
SI  1= AKSD# *2*  < GAMD* *2-7 **Z ) 

S22  =  AK3D**2* <  AKAP**2-Z**2 ■  '3  .  M  GAMD*Z ' **2 
DELS=S11*S22-<ANSD*Z*<  GAMD  **2-2.  >>.**.? 

F  =  CSQRT  <  <Z«AKFD/AKSD;**2-1.  . 

GO  TO  (  10,20* 20*30,30'  ID 
DSP  =  AI#DELA#DELS-Z**3*AN3D#  •  RHO/2  .  '  *  ■  DEL**  ‘Ell  rDE'_'-*Al  1 
GO  TO  40 

DSP=AI*DELA-Z*#3*AKSD*RH0*A1  1/ <F«FF  ) 

GO  TO  40 

DSP  =  AI*DELS-Z  K*3#AKSD*RH0*S1 1 /<F#FF  ' 

RETURN 
END 


I 
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FUNCTION  DSF‘2  (  Z  ) 

COMPLEX  DSF'2.  AI  .  All  .  A22.S11  »  S22  *  AKD  *  SF'R  *  GF'R  *  AKAP  *  D 1  *D2* 
COMPLEX  DELA*DELS.Z»F*CMPLX.C3GRT*CC0S*CSIN 
COMPLEX  CCQ.CSG* CSS. CCS 

COMMON  RHO.  AKSD . GAMP . GAMD *  AKAR.AKAI  *  AKFD *  ID*  IDoF 
IBD=ID-5 
F  F  =  2  . 

GO  TO  <  7  *  6  >  5  *  6  *  5  )  IBD 

5  FF=  1  . 

GO  TO  7 

6  F  F  -  2  . 

7  CONTINUE 
AI=CMF'LX<0.  .1  .  ) 

AKD= AKSD/Z 

IF  (REAL<AK0#*2-AKSD**2> >  11.12*12 

12  SPR=-AI*CSQRT< AKB**2-AKSB**2 ) 

GO  TO  13 

1 1  SF'R  =  CSGRT<  AKSD**2-AKDf  *2) 

13  IF<REAL(1.-<Z/GAMD)*#2))  21*22*22 

21  QF'R  =  CSQRT  (  (  Z/GAMO  )  **2-1  .  )  *  A  K  D 

GO  TO  23 

QF'R  =  -AI*CSQRT(1  ♦ - < Z/6AMP ) *#? > *AKD 
A 1 1  =  (  A  K  D  *  *  2  -  S  F'  R  *  *  2  )  *  C  3 1 N  \  Q  P  R  ) 

A 2 2  =  -  (  AKD**2-SPR**2  )  *CCOS  <  SF'R  > 

DELA-A1  l*A22-4  ♦  *AKD**2*GPR*SPR*C3  f  N  *'  3 s- R  •  KCCl  -v  (  OPR  ) 

SI 1  =  < AKD**2-SPR**2)*CCGS •OPR) 

S22  =  -(AKH**2-SF’R**2)*CSIN(SFf\) 

DELS*S1 1  *S22-4  .  *AKD**2*GPR*SPR  *£:S  IN  <  QPR  >  *CC(-S  •  SF'R  ) 

F  =  CSGRT ( <Z*AKFD/AKSD)**2-1  ,  ) 

D1=DELA*DELS*AI 

D2«RH0*QPR*AKSD**4*DELS*CC0S  <  GF’R  )  *CCOS  <  S='R  )  /  AKD*F  > 

D3  =  RH0*GF'R*AKSD**4*DELA*CSIN<  QPR  )  *  C  S  I  N  <  S  P  R  »  ,  <  A  N  0  *  F  ) 

GO  TO  <15*25*25*35*35)  I  DC 
DSF'2  =  D 1  +  (  D3-D2 )  /FF 
GO  TO  45 
DSF'2  =  D  1 -D2/FF 
GO  TO  45 
D  3  P  2  =  D 1  +  D  3  /  F  F 
CONTINUE 
RETURN 
END 
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APPENDIX  E 


EXAMPLES 


Items  1  and  2  show  two  examples  where  the  running  of  the  program  took 
place  without  complications.  In  Item  3  the  program  fails  to  converge, 
conceivably  due  to  the  close  proximity  of  another  root  of  the  dispersion 
relation,  originating  in  a  symmetric  real  root.  This  example  shows  how 
increasing  ANR,  which  is  the  number  of  steps  to  reach  the  full  value  of 
the  density  of  the  fluid,  does  result  in  convergence  of  the  program. 


RUN  R 0 8 1 9 

1  ENTER  POI  .COCS  .  RHM  I  .3028.  .49332.  .1280* 

ENTER  A R A R  #  ANA  I »  DELZ  I  .9278.0.0,0.01 
ENTER  ANSD.ZNOT.*  OF  IT. :4. 05135* .9003,200. 

ENTER  ANR.  TOL.  FSTEF :  1., 0.01,1. 

ENTER  SID!  2- 
DEL-  0 . 33621 E-22 

RHO.REGA.AMGA.DELH*  0.128060  0.900309  0.33512E-02  01AB10E-07 

THICK  PLATE.  A. ONLY,  FLUID  ONE  SIDE 
UALUES  OF  ANSD.ZNOT.  ARE  9.051350  0.900300 

REAL  AND  IM.  PARTS  OF  GAMMA-ZNOT  AND  ERROR  ARE 
0 . 90003E-05  0.33512E-02  0.16810E-07 

DIMENSIONLESS  ATTEN. FACTOR  =  0.16750E--01 

ENTER  RETURN  OPTION  NUMBER  3 

-  ENTER  ANSD.ZNOT,*  OF  IT . 15.9923, .9192,200. 

ENTER  ANR,  TOL,  FSTEP!  1., 0.01.1. 

ENTER  BID!  2. 

DEL-  0. 23992E-02 

RHO.REGA.AMGA.DELH*  0.128060  0.919153  0.23979E-02  0.11971E-06 

THICK  PLATE.  A. ONLY.  FLUID  ONE  SIDE 
UALUES  OF  AKSD.ZNOT,  ARE  5.992300  0.919200 

REAL  AND  Ih.  PARTS  OF  GAMMA-ZNOT  AND  ERROR  ARE 
-0 . 97207E-09  0.23979E-02  0.11971E-06 

DIMENSIONLESS  ATTEN. FACTOR  =  0.17051E-01 

ENTER  RETURN  OPTION  NUMBER  3 

3  ENTER  ANSD.ZNOT,*  OF  I T . ! 7 . 8 1 795 . . 9 1 97 . 200 . 

ENTER  ANR,  TOL.  FSTEP :  1., 0.01.1. 

ENTER  BID:  2. 

DEL  =  0. 18673E-02 
EXIT  LEFT 

ENTER  RETURN  OPTION  NUMB'ER  9 

ENTER  ANR.  TOL.  FSTEP:  3. ,0.01,1. 

ENTER  BID:  2. 

DEL  -  0 . 62295E-03 

RHO.REGA.AMGA.DELH*  0.092687  0.919709  0.61999E-03  0.31122E-07 
RHO.REGA.AMGA.DELH*  0 . 0853  73  0.919708  0.12900E-02  0.3U22E-07 
EXIT  RIGHT 

ENTER  RETURN  OPTION  NUMBER  9 
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enter  anr. 

ENTER  BID:  2. 

DEL  *  0.37347E-03 
RHQ  , REGA »  AMGA  » DELH* 
RHO.REGA.ftNGA.DELH* 

rho.rega. amga.delh* 

RHO.REGfl. ANGfl. DELH* 
RHO.REGA. AMGA.DELH* 
THICK  RLATE.  A. ONLY. 
VALUES  OE  AKSD.ZNOT. 


TOL.  ESTEP!  3. .0.01,1. 


0.023612  0.919709  0.37199E-03 

0.031224  0.919709  0.74399E-03 

0.076836  0.919708  0.11160E-02 

0.102448  0.919707  0.14880E-02 

0.128060  0.919705  0.18600E-02 

FLUID  ONE  SIDE 
ARE  7.317450  0.919700 

REAL  AND  IN.  PARTS  OF  GANMA-ZNOT  AND  ERROR  ARE 
0.50664E-05  0.19600E-02  0.18673E-07 

DIMENSIONLESS  ATTEN. FACTOR  *  0.17190E-01 

ENTER  RETURN  OPTION  NUMBER  0 


0. 18673E-07 
0. 18673E-07 
0. 18673E-07 
0. 1B673E-07 
0. 18673E-07 
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The  next  example  shows  that  in  some  cases  increase  of  ANR  or  increas¬ 
ing  FSTEP  can  make  the  program  converge.  It  is  not  clear  whether  or  not 
there  are  cases  where  one  of  these  two  variables  would  induce  convergence 
and  not  the  other. 


RUN  R 0 8 1 9 

ENTER  POI .CQCS.RHM:  . 3028 .. 4432 ,  . 1 2806 
ENTER  AKAR,  AK  A I  ,  DEL2:  . 9278 , 0 . 0 , 0 . 0 1 
ENTER  ANSD  » ZNOT  « *  OF  I T . : 7 . 5249 , . 9290 » 200 . 
ENTER  ANR.  T0L.  FSTEP*.  3. .0.01.1. 

ENTER  BIH:  6. 


DEL=  0 . 54452E-03 

RHO.REGA, AMGA»DELH=  0.042687  0.927811  0.64525E-03 

RHO.REGA. AMGA»DELH=  0.085373  0.927847  0.26817E-03 

RHO.REGA, AMGA,D£LH=  0.128060  0.927850  0.17425E-03 

EXACT.  A.  AND  S..  FLUID  ONE  SIDE 
VALUES  OF  AKSD.ZNOT.  ARE  7.524900  0.929000 

REAL  AND  IM.  PARTS  OF  GAMMA-ZNOT  AND  ERROR  ARE 
-0 . 1 1500E-02  0 . 1 7425E-03  0.27226E-06 

DIMENSIONLESS  ATTEN. FACTOR  =  0.15230E-02 

ENTER  RETURN  OPTION  NUMBER  5 


0 . 27226E-05 
0 . 27226E-06 
0 . 27226E-06 


ENTER  BID!  9. 

DEL  =  0 . 58163E-03 
EXIT  LEFT 

ENTER  RETURN  OPTION  NUMBER  4 


ENTER  ANR,  TOL,  FSTEP :  10. ,0.01,1. 

ENTER  BID:  9. 

DEL=  0 . 1 7449E-03 

RHO,REGA,AMGA,DELH= 

0.012806 

0.929164 

0.47871E-03 

0 . 87244E-07 

RHO.REGA, AMGA,DELH= 

0.025612 

0.929167 

0 . 95742E-03 

0 . 87244E-07 

RHO,REGA,AMGA,DELH= 

0.038418 

0.929170 

0. 14361E-02 

0 . 87244E-07 

RHO.REGA, AMGA,DELH= 

0.051224 

0.929175 

0. 19148E-02 

0.87244E-07 

RHO.REGA, AMGA,DELH= 

0.064030 

0 . 929182 

0 . 2393BE-02 

0 . 87244E-06 

RHO.REGA, AMG A, DELH= 

0.076836 

0.929190 

0 . 28728E-02 

0. 87244E-06 

RHO.REGA. AMGA,0ELH= 

0.089642 

0.929200 

0 . 335 1 7E-02 

0 . 87244E-06 

RHO.REGA, AMGA,DELH= 

0.102448 

0.929212 

0 . 38290E-02 

0. 87244E-06 

RHO.REGA. AMG A, DELH= 

0.115254 

0.929224 

0. 43079E-02 

0 . 87244E-06 

RHO.REGA, AMGA,DELH= 

0.128060 

0.929237 

0 . 47869E-02 

0 . 87244E-06 

EXACT,  S.  ONLY,  FLUID  ONE  SIDE 

VALUES  OF  AKSD.ZNOT 

,  ARE  7.5 

24900  0.929000 

REAL  AND  IM.  PARTS 

OF  GAMMA-ZNOT  AND  ERROR 

ARE 

0 . 237 1 1 E-03  0 

. 47869E-02 

0 . 87244E- 

06 

DIMENSIONLESS  ATTEN 

.FACTOR  = 

0. 41715E-01 

ENTER  RETURN  OPTION 

NUMBER  4 

ENTER  ANR,  TOL,  FSTEP:  3.  ,0.01 

,3. 

ENTER  BID:  9. 

DEL  =  0 . 1 7449E-02 

RHO.REGA, AMGA.DELH* 

0.042687 

0.929172 

0. 15958E-02 

0 . 87244E-07 

RHO.REGA, AMGA»DELH= 

0.085373 

0 . 929196 

0.31915E-02 

0 . 87244E-06 

RHO.REGA, AMGA,DELH= 

0.128060 

0.929237 

0 . 47872E-02 

0 . 87244E-06 

EXACT,  S.  ONLY.  FLUID  ONE  SIDE 

VALUES  OF  AKSD.ZNOT 

,  ARE  7.5 

24900  0.929000 

REAL  AND  IM.  PARTS 

OF  GAMMA-ZNOT  AND  ERROR 

ARE 

0.23729E-03  0 

. 47872E-02 

0 . B7244E- 

06 

DIMENSIONLESS  ATTEN 

.FACTOR  = 

0.4171 7E-01 

ENTER  RETURN  OPTION  NUMBER  0 
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